Q1¶

$1$. A random variable $X$ is said to follow gamma$(\alpha,\beta)$ distribution.If it's probability density function given by $$\begin{equation} f(x) = \begin{cases} \frac{e^{-\frac{x}{\beta}}x^{\alpha-1}}{\beta^{\alpha}\Gamma(\alpha)}, 0 < x< \infty \\ 0~~~~~~ ,otherwise \end{cases} \end{equation}$$.

(a) . Draw $f(x)$ for different values of $\alpha$ and $\beta$. Create two plots. In one plot fix $\alpha = 2$ and vary $\beta = 2,4,6.5$. In the other plot, fix $\beta = 3$ and vary $\alpha = 4,6,7.5$. Each plot will contain three functions. Use different functions.Add appropiate legends to the plots.

(b). It is known that if $X_{1},X_{2},\cdots,X_{m} \sim$ gamma$(1,\beta)$, then $\sum_{i = 1}^m X_{i} \sim \mathcal{gamma}(m,\beta)$. Verify the statement using simulation. Use $m = 10$ and $\beta = 3$. Consider $1000$ replication to visualize the distribution of $\sum_{i =1}^m X_{i}$.

In [28]:
using Distributions,Plots,SpecialFunctions,QuadGK;
In [2]:
using LaTeXStrings;
In [29]:
function f(y,a,b)
    (exp(-y/b)*(y^(a-1)))/((b^a)*gamma(a))
end
Out[29]:
f (generic function with 3 methods)
In [30]:
a = 2
b_vals= [2,4,6.5];
y = 0:0.01:15;
p = plot()
for b  in b_vals
    l = [f.(y,a,b)]
    p =  plot!(y,l,label = "Gamma($(a),$(b))",fg_legend = false)

end
display(p)
In [35]:
b = 3
a_vals=  [4,6,7.5]
y = 0:0.01:15;
p1 = plot()
for a  in a_vals
    l = [f.(y,a,b)]
    p1=  plot!(y,l,label = "Gamma($(a),$(b))",legend = :topleft,fg_legend = false, grid = false)
end
display(p1)
In [9]:
In [ ]:

Q 1 b¶

In [21]:
x1 = []
m = 10
beta = 3
rep = 1000
for i in 1:rep
    push!(x1,rand(Gamma(m,beta)))
end
In [25]:
# beta density function with (m,beta) as parameters
function f(x)
    (exp(-x/beta)*x^(m-1))/(gamma(m)*(beta^m))
end
Out[25]:
f (generic function with 2 methods)
In [ ]:
In [31]:
histogram(x,normalize = true,label = "histogram of Xi",c = 4,fg_legend = false)
plot!(f,10,60,label = "gamma(10,3)",c=  6, linewidth = 4,style = :dash)
Out[31]:

Q2¶

$2$. Consider the following function $\begin{equation} g(x) = \begin{cases} \frac{1}{\sqrt{2x}}e^{-\frac{x^2}{4}} , 0 < x< \infty \\ 0 ~~~~~~~~, otherwise \end{cases} \end{equation}$

(a) Using Monte Carlo method compute the integral of the function on the interval $[a,b] = [2,5]$ based on $n = 1000$ samples. Call it $I_{n}$.

(b) Using inbuilt function to compute the integral, call it $I$.

(c) Print the absolute difference between $I$ and $I_{n}$.

(d) Try simulate, show that $I_{n}$ converges to $I$ as $n \to \infty$.

We want to compute the integral $\int_{2}^{5} g(x) dx$.Take $a = 2,b = 5$.Then the integral can be written as $(b-a)\int_{a}^{b}\frac{1}{b-a}\frac{1}{\sqrt{2x}}e^{-\frac{x^2}{4}}dx = E_{f}[g(X)] ,f\sim \mathcal{U}(a,b)$.

So we generate random samples $x_{i}$ from $\mathcal{Uniform(a,b))}$, then calculate $g(x_{i})$ then take the mean and multiply it by $(b-a)$. We do this process for $n = 1,2,\cdots, 1000$.

In [10]:
function g(x)
    (1/sqrt(2*x))*(exp(-x^2/4))*(x>0)
end
Out[10]:
g (generic function with 1 method)
In [13]:
a = 2; b = 5;
I_n = []
n = 1000;

for i in 1:n
    x = rand(Uniform(2,5),i)
    push!(I_n,(b-a)*mean(g.(x)))
end
In [ ]:
In [14]:
I_n
Out[14]:
1000-element Vector{Any}:
 0.007956155850635398
 0.012879270597985045
 0.16372531948611774
 0.10490097549397108
 0.018840875047610554
 0.18014317027367047
 0.12268589027643673
 0.2560520985104784
 0.11287763607330115
 0.10772582282532644
 0.1291242512115906
 0.1421804375174541
 0.10591406591090198
 ⋮
 0.1159685271251723
 0.12555332242330125
 0.12746130845837733
 0.13371915904913068
 0.12781254015601573
 0.12663111881770747
 0.1242725332912191
 0.12648038568286088
 0.12775120893860106
 0.12413690244898444
 0.12225380072351498
 0.12431813357814087
In [15]:
mean(I_n)
Out[15]:
0.12245044474453692
In [16]:
exact_integral = quadgk(g,a,b)[1]
Out[16]:
0.12290701094708961
In [23]:
plot(I_n,label = "approx integral")
hline!([exact_integral],label = "exact integral")
Out[23]:
In [21]:
abs_diff = []
for j in 1:length(I_n)
    push!(abs_diff,I_n .- exact_integral)
end
length(abs_diff)
Out[21]:
1000
In [27]:
plot([abs_diff],legend = false)

hline!([0],label =:none)
Out[27]:
In [ ]: